
    #################################################################
    ### Figure 4 (Senate) - R code

    T1 <- 92
    T2 <- 114
    T <- length(T1:T2)  # number of congresses
    year.odd <- seq(from=1971, to=2016, by=2)
    

    ### Import the individual-level data

    Micro_Data1 <- read.csv("Figure4_Table2_Data_Senate.csv", header=T)
    dim(Micro_Data1)
    names(Micro_Data1)
    
    Micro_Data <- na.omit(Micro_Data1)
    dim(Micro_Data)
    names(Micro_Data)
    
 
    ### Get the Congress-level data
    Macro_Data <- read.csv("Table_2_3_Congress_level_data.csv", header=T)
    names(Macro_Data)

    library(mvtnorm)
    library(R2WinBUGS)
    
    ### Setup data
    Y  <- Micro_Data$EP.extreme 
     
    Party <- ifelse(Micro_Data$party==100, 1, 2)
    GOP <- ifelse(Micro_Data$party==200, 1, 0)
    Southern.State <- c(40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 51, 54) # VA, AL, AR (Arkansas), FL, GA, LA, MS, NC, SC, TX, plus KY and TN

    X1 <- Micro_Data$Energy_cmt
    X2 <- Micro_Data$Interiror_cmt
    X3.1 <- Micro_Data$Environment.extreme
    X3  <- as.vector((X3.1 - mean(X3.1, na.rm=T))/(sd(X3.1, na.rm=T)))

    Divided <- ifelse(Macro_Data$President_party==Macro_Data$Senate_Majority, 0, 1)
    Majority_Margin <- Macro_Data$Senate_Majority_Seats/Macro_Data$Senate_Total_Seats
    gao_reports <- Macro_Data$gao_reports
    Oil <- log(Macro_Data$oil_real) 
    Warming <- Macro_Data$global_warming
    
    # Standardize the continuous variables
    Majority_Margin.std  <- as.vector((Majority_Margin - mean(Majority_Margin, na.rm=T))/(sd(Majority_Margin, na.rm=T)))
    gao_reports.std  <- as.vector((gao_reports - mean(gao_reports, na.rm=T))/(sd(gao_reports, na.rm=T)))
    Oil.std  <- as.vector((Oil - mean(Oil, na.rm=T))/(sd(Oil, na.rm=T)))
    Warming.std  <- as.vector((Warming - mean(Warming, na.rm=T))/(sd(Warming, na.rm=T)))
    Warming.std.lag  <- c(Warming.std[1], Warming.std[-length(Warming.std)]) 
    
    Z1 <- Divided             
    Z2 <- Majority_Margin.std 
    Z3 <- Oil    
    Z4 <- Warming.std.lag
    Z5 <- gao_reports.std
  
    congress <- Micro_Data$time

    # create senator indicators
    ICPSR_ID <- as.vector(Micro_Data$ICPSR_ID)
    uniq.icpsr_id <- unique(ICPSR_ID)
    n.senator <- length(uniq.icpsr_id)
    senator <- rep(NA, n.senator)
        for (i in 1:n.senator){
            senator[ICPSR_ID == uniq.icpsr_id[i]] = i
        }

    K <- n.senator
    N <- length(Y)
    J <- length(Z1)
    
    n.gamma <- 6 # number of macro-level coefficients
   
    data <- list("N", "J", "K", "n.gamma", "Y", 
                 "X3", "congress", "senator", 
                 "Z1","Z2","Z3","Z4","Z5")

    # setup initial values
    init1 <- list(beta1=rnorm(J), sigma.y=runif(1),
                  gamma=rnorm(n.gamma), alpha1=rnorm(K), 
                  alpha=rnorm(J), sigma.alpha=runif(1))
    init2 <- list(beta1=rnorm(J), sigma.y=runif(1),
                  gamma=rnorm(n.gamma),alpha1=rnorm(K),
                  alpha=rnorm(J), sigma.alpha=runif(1))
    init3 <- list(beta1=rnorm(J), sigma.y=runif(1),
                  gamma=rnorm(n.gamma),alpha1=rnorm(K),
                  alpha=rnorm(J), sigma.alpha=runif(1))

    inits <- list(init1, init2, init3)
    
    # setup parameters
    parameters = c("beta1", "sigma.y", "alpha1", 
                    "gamma", "alpha", "sigma.alpha")
    
    Model.fit = bugs(data, inits, parameters, "figure4_senate.bug",
                     n.chains=3, n.thin=10, n.burnin=2000, n.iter=12000)
    
    ### 4. Analyze the draws from the Posterior distribution ############################################

    library(coda)    
    chain1 <- read.coda(output.file="coda1.txt", index.file="codaIndex.txt", quiet=T)
    chain2 <- read.coda(output.file="coda2.txt", index.file="codaIndex.txt", quiet=T)
    chain3 <- read.coda(output.file="coda3.txt", index.file="codaIndex.txt", quiet=T)

    MCMC <- rbind(chain1, chain2, chain3) # Choose the chain that puts Democrats on the left

    iter <- nrow(MCMC)

    alpha.finder  <- seq(from=1, to=T, by=1)
    alpha1.finder  <- seq(from=(T+1), to=(T+K), by=1)
    beta1.finder <- (T+K+1):(T+K+T)
    gamma.finder  <- seq(from=(T+K+T+1+1), to=(T+K+T+n.gamma+1), by=1)

    n.beta1 <- length(beta1.finder)

    alpha.mc  <- MCMC[,alpha.finder]
    beta1.mc   <- MCMC[,beta1.finder]
    gamma.mc  <- MCMC[,gamma.finder]

    alpha.label <- paste(T1:T2, "th Congress", sep="")
    beta1.label <- year.odd
    
    gamma.label <- c("Intercept                                       ", 
                     "Divided Government",            
                     "Majority Party Seat Margin",
                     "Oil Price (log)", 
                     "Global Warming",
                     "GAO Reports"
                     )

    ### (1) Histograms of marginal posteriors of selected parameters
    ################################################################

    postscript("density_beta1_trace.eps", paper='letter', horizontal=F)
    par(mfrow=c(6,2), mar=c(2, 5, 3, 3)) # 'mar -> c(bottom, left, top, right)'
    for(i in 1:6){
    plot(density(beta1.mc[,i]), main="", xlab="")
    title(main=paste("Marginal Posterior:", beta1.label[i]), font.main=1)
    plot(beta1.mc[,i], type="l", xlab="iteration", ylab="")
    title(main="Traceplot", font.main=1)
    }
    dev.off()

   
    ### Credible Intervals
    beta1.CI  <- round(apply(beta1.mc,  2, FUN=quantile, probs=c(0.025, 0.975, 0.5, 0.05, 0.95, 0.005, 0.995)), 3)
    gamma.CI  <- round(apply(gamma.mc,  2, FUN=quantile, probs=c(0.025, 0.975, 0.5, 0.05, 0.95, 0.005, 0.995)), 3)
    alpha.CI  <- round(apply(alpha.mc,  2, FUN=quantile, probs=c(0.025, 0.975, 0.5, 0.05, 0.95, 0.005, 0.995)), 3)

   
    plan <- rbind(c(1, 1), c(2, 2))
    print(plan)

    postscript("figure4_senate.eps", height=7, width=5, horizontal=F)
    layout(plan)
    par(mar = c(0,0,0,0), oma = c(3,4,0,0), mgp=c(1.5, 0, 0), xpd=NA, tck=-0.007)  #c(bottom, left, top, right)

    plot(1:n.beta1, 1:n.beta1, type='n', col='black', ylim=c(-1, 1.5), 
            xlab="", ylab="", xaxt="n")
    text(    1:n.beta1, beta1.CI[5,1:n.beta1]-0.5, labels=beta1.label,  cex=.8, srt=90) # use srt for text rotation 90=degree
    points(  1:n.beta1, beta1.CI[3,1:n.beta1], type='p', pch=19, col='blue') # Plot the Medians
    segments(1:n.beta1, beta1.CI[1,1:n.beta1],  
             1:n.beta1, beta1.CI[2,1:n.beta1], col='red', lwd=2) # Plot 95% CI 
    segments(1:n.beta1, beta1.CI[4,1:n.beta1],  
             1:n.beta1, beta1.CI[5,1:n.beta1], col='black', lwd=3) # Plot 90% CI
    title(main="",  font.main=1)
    points(  0:n.beta1, rep(0, n.beta1+1), type='l', lty=2) # use this instead of abline 

    plot(1:n.gamma, 1:n.gamma, type='n', col='black', ylim=c(-1, 1.5), 
           xlab="", ylab="", xaxt="n") # , yaxt="n"
    text(    1:9 - 0.15, gamma.CI[5,], labels=gamma.label,  cex=1, srt=90) # use srt for text rotation 90=degree
    points(  1:n.gamma, gamma.CI[3,1:n.gamma], type='p', pch=19, col='blue') # Plot the Medians
    segments(1:n.gamma, gamma.CI[1,1:n.gamma], 
             1:n.gamma, gamma.CI[2,1:n.gamma], col='red', lwd=2) # Plot CI
    segments(1:n.gamma, gamma.CI[4,1:n.gamma], 
             1:n.gamma, gamma.CI[5,1:n.gamma], col='black', lwd=3) # Plot CI
    #title(main=paste("Congress-Level Coefficients"),  font.main=1)
    points(  1:n.gamma, rep(0, n.gamma), type='l', lty=2) # use this instead of abline 
    dev.off()
    
